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Abstract 

Soliton rate equations are based on non-Kolmogorovian models of proba- 
bility and naturally include autocatalytic processes. The formalism is not 
widely known but has great unexplored potential for applications to systems 
interacting with environments. Beginning with links of contextuality to non- 
Kolmogorovity we introduce the general formalism of soliton rate equations 
and work out explicit examples of subsystems interacting with environments. 
Of particular interest is the case of a soliton autocatalytic rate equation cou- 
pled to a linear conservative environment, a formal way of expressing sea- 
sonal changes. Depending on strength of the system-environment coupling 
we observe phenomena analogous to hibernation or even complete blocking 
of decay of a population. 
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1. Introduction 



It is quite typical that soliton equations discovered in one domain of 
science turn out to have apphcations in completely different fields. The so- 
called nonlinear Schrodinger equation is a perfect example. Its applications 



range from surface waves on deep water (Zakharov, 1968 ) to nonlinear effects 



in DNA (Bishop et al. 1980) and fiber optics (Kibler et al. 2010). The 



name comes from mathematical similarity to the basic equation of quantum 
mechanics. However, this is only a similarity. The true Schrodinger equation 
is always linear. 

In this paper we concentrate on another soliton system, Darboux-integrable 



von Neumann equations (Leble & Czachor, 1998; Ustinov et al. 2001 Cieslihsd 



Czachor & Ustinov, 2003), whose origin goes back to studies on generaliza- 



tions of quantum mechanics (Czachor, 1997), but which has never gained 



great popularity among physicists. The reasons are similar to those from 
our previous example: von Neumann equations occurring in quantum me- 
chanics are always linear, but Darboux-integrable von Neumann dynamics 
can be nonlinear. It naturally includes catalytic and autocatalytic processes 



(Aerts & Czachor, 2006 Aerts et al. , 2006), but the probability calculus be 



hind it is contextual and hence non-Kolmogorovian. So, is it possible that 
the formalism discovered in attempts of generalizing quantum mechanics has 
unexpected applications in a different domain? If so, what kind of criteria 
should one use to identify these new applications? The hints come from 
autocatalysis and probability. 

Contextuality plays in probability a role similar to that of curvature in ge- 
ometry. In geometry curvature measures non-commutativity of translations. 
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In probability contextuality measures non-commutativity of questions. In 
geometry curvature requires local charts: taken together they form a global 
object, an atlas, covering the entire curved space. The charts are pairwise 
compatible so that one can explore the whole curved object by thumbing the 
atlas. The whole global structure is called a manifold. In probability, ques- 
tions asked earlier define contexts for those that are asked later. Typically, 
a triple of questions A ^ B ^ C has to be described in two steps, A ^ B 
and then B C, and each step involves a different probability space. The 



global structure is not a probability space in the sense of Kolmogorov (1956 ). 



Rather, it is an object of a manifold type (Gudder, 1984). 

Non-Kolmogorovian probabilities are sometimes termed quantum proba- 
bilities, but this name is misleading and its origins are historic and not on- 
tological. Non-Kolmogorovian structures are as ubiquitous as non-Euclidean 



geometries (Khrennikov, 2010). Non-Euclidean geometries were discovered 
before the advent of general relativity, otherwise one would speak of "gravi- 
tational geometry" . Non-Kolmogorovian probability was discovered in quan- 
tum mechanics, but after some 30 years of studies in logical foundations of 
quantum mechanics it has become clear that non-Kolmogorovity has noth- 



ing to do with microphysics (Accardi & Fedulo, 1982; Aerts, 1986 Pitowsky 



1989). It is thus striking that the search for fundamental laws of ecology 



has led ecologists to probabilistic structures of a propensity type (Ulanowicz 



1999, 2009), which are known to be conceptually close to probability models 



of quantum mechanics (Popper 1982). 

The goal of the present paper is to introduce formal mathematical struc- 
tures that are needed in discussion of nonlinear models based on non-Kolmogorovian 



3 



probability. In Section 2 we begin our analysis with a simple illustration of 
non-Kolmogorovity: A cyclic competition. In Section 3 we discuss a model 
of non-Kolmogorovian probability (based on vectors), and further generalize 
it in Section 4 to the density-matrix formalism. In Section 5 the density- 
matrix formalism is used to reformulate (possibly nonlinear) rate equations 
in a Lax-von Neumann form. In order to develop some intuitions we first 
concentrate on a linear example. In Section 6 we introduce the notion of 
a hierarchy of coupled environments. The first examples of nonlinear rate 
equations and their Lax-von Neumann representation occur in Section 7. Sec- 
tion 8 introduces the main subject of this paper: soliton rate equations. An 
equation describing a linear environment coupled to a nonlinear subsystem is 
explicitly analyzed. Choosing various values of parameters we plot evolution 
of populations for shorter and longer time scales and alternative couplings 
between systems and their environment. All these examples are based on 
exact solutions of the associated system of nonlinear rate equations. Section 
9 is devoted to the specific case of a periodically changing environment. Fi- 
nally, in Section 10, we give a preliminary analysis of soliton systems with 
dissipation and the role of environment in possible blocking or hibernation of 
dissipative decay of populations. Section 11 is the first step that goes beyond 
the soliton dynamics. Here we introduce a new generalization of replicator 
equations applicable to games where players try to modify the rules of the 
game during its course. The equation is of a von Neumann form and seems 
to have applications to real-life evolutionary games, a subject that will be 



discussed in more detail in a future work, cf. (Aerts et al. , 2013). 
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2. Cyclic competition is non-Kolmogorovian 

A simple form of cyclic competition occurs in the rock-paper-scissors 
game: Rock R destroys scissors S, scissors destroy paper P, paper destroys 
rock. In a kinetic form the game corresponds to 

R + S — > R + decay products, (1) 
S + P —7- S + decay products, (2) 
P + R —7- P + decay products. (3) 

Logical and probabilistic version of the game is illustrated by Fig. 1. There 
are two players and three random variables P, R with binary values 0, 1. 
Players choose which random variables will be measured, and then appropri- 
ate pairs of binary values occur with probabilities p{S = 0, R = 1) = p{P = 
0, S = 1) = p{R = 0, P = 1) = 1. The remaining probabilities vanish. 

In spite of trivial statistics the problem is formally quite subtle, a fact 
mentioned in the context of "non-quantum" probability already by [Vorobev 



(1962). Vorobev's ideas led mathematicians to the notion of a contextual 



marginal problem (Fritz & Chaves 2012), i.e. the question if and when a 



collection of probabilities can be regarded as a result of computing marginal 
probabilities for pairs from joint probabilities of triples. If triple joint prob- 
abilities do not exist, one has to resort to local probability spaces, and one 
arrives at a manifold-like non-Kolmogorovian structure. 

In fact, we can model the game on three independent probability spaces 
corresponding to the three pairs of random variables shown in Fig. la. Each 
such pair can be realized in an experiment. The resulting data can be ma- 
nipulated in standard ways, but one has to be cautious not to mix data from 



different alternative experiments. To give an example, the average of S can 

be computed by means of 
1 

S=J2 MS = a,P = b) =p{S = 1,P = 0) = 1. 

a,b=0 

If one computes the average as follows, 

1 

S = ap{S = a, R = b) = 0, 

a,b=0 

a contradictory result is obtained. This does not mean that = 1. There is no 
contradiction if one takes into account that logically it is not allowed to treat 
S in the S-P experiment as the same 5* as in the S-R one. R and P form 
different contexts for S. In quantum mechanical terminology one can say 
that the two players constitute a non-local system. One can also reformulate 
the RPS game in a way that explicitly turns it into a variant of Einstein- 



Podolsky-Rosen experiment (Aerts et al. , 2011) where Bell's inequality (Bell 



1964) is violated. 



But then we arrive at another difficulty. The rock-paper-scissors game 
is a canonical example of cyclic competition. It has its analogues in pop- 
ulation dynamics of plankton, lizards or bacteria. Standard approaches to 
kinetic dynamics imply that elementary processes ([l])-(|3]) lead to nonlinear 
rate equations of the form 

[S] = -ks[S][R] + . . . , (4) 
[R] = -kn[R][P] + ..., (5) 
[P] = -kp[P][S] + .... (6) 

The question is if [S] in Q is the same [S] as the one in (|6])? If so, why? 
Clearly, Q represents a process where S interacts with R, while ^ occurs 



due to interactions between P and S. But we have just shown that a naive 
mixing of two the two Ss leads to contradictions. 

So the difficuhy is that we have to combine simultaneously two apparently 
contradictory aspects. On the one hand, in the RPS game the other player 
chooses either P or R, but P and R cannot be chosen simultaneously. The 
two contexts for S cannot occur at the same instant of time. On the other 
hand, however, [S] in both equations is a time dependent variable [5'](t), 
with the same values of t in Q and ([6]). The problem is typical for quan- 
tum mechanical evolutions. Its solution is known since 1920s, when it was 
understood how to probabilistically interpret solutions of the Schrodinger 
equation. 

3. Non-Kolmogorovian probability 

Kolmogorovian model of probability leads to conceptual difficulties when- 
ever order of "questions" (i.e. context) is not irrelevant. The Kolmogoro- 
vian algorithm for constructing conditional and joint probabilities is formally 
based on projections on subsets of events, Xa^ = A C X. In practical ap- 
plications xa is often represented by characteristic functions 



If yU is a probability measure on the set of events X, n{X) = 1, then proba- 
bility of finding x & A G X is defined as 





p{A) 



KxaX). 



(8) 
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In case of joint probabilities 

p{A nB) = nixAnsX) = ij,{xaXbX) 

= fi{xBXAX)=p{BnA). (9) 

Conditional probabilities are defined via the Bayes rule 

p{AnB) = p{A)p{B\A), (10) 
pjAnB) ii{xaXbX) 

Based on set-theoretic properties one can derive a number of inequalities that 
have to be satisfied by p{A). For example, < p{A) < 1, p{A H B) < p{A), 
or 

p{AnBr\C)=p{Ar\Cr\B) <p{AnB). (12) 



One of the simplest problems with (12) occurs in quantum mechanics with 



the so-called Malus law stating that joint probability of photon's passage 
through two polarizers A and B is 

p{A n B) = p{A)p{B\A) = - cos'^ (j) (13) 

where is the angle between the polarizers and p{A) = 1/2 is the proba- 
bility of transmission through the first polarizer. Obviously, if A and B are 
perpendicular then p{A (1 B) = since p{B\A) = 0. But for polarizers tilted 
by 45 degrees one finds p{B\A) = 1/2. In consequence, if one takes a third 
polarizer C tilted by 45 degrees with respect to mutually perpendicular A 
and B, then 

p{AnCnB) = p{A)p{C\A)p{B\C) 
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p{AnBn C) 



p{A)p{B\A)p{C\B) 




(15) 



p{AnB) = 0, 



(16) 



leading to the counterintuitive inequality 



p{AnCn B) 



1 



p{AC\B). 



(17) 
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> 



By ^4 n C n S we denote an event where the photon first crosses A, then C, 
and finally B. From a logical point of view the interpretation makes sense. 
Indeed, if the source of light is to the left of A and the detector is to the right 
of B, then the act of detection means that the particle had to be transmitted 
through all the polarizers: It had to pass through A, and C, and B. It follows 
that the detection event can be regarded bs AnC HB. 

Symbolically, for photons and polarizers A ^ C ^ B may be allowed, 
whereas the direct process A ^ B may be forbidden. This is the essence of 
the problem and not the fact that we speak of quantum particles. Analogous 
diagrams are typical of chemistry, biology, psychology, cognitive science... 

In such theories, ii p{Ar\C r\B) ^ p{Ar\Br\C) then one cannot assume 
that 



One needs a different mathematical model. In some (but not all) generalized 
models of probability one begins with the basic property of "questions" (or 
"propositions" in the logic parlance) , 



p{Ar\Br\C)^ f^iXAnBncX) = fxixAXBXcX). 



(18) 



XAnA = Xa = XaXa, 



(19) 
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representing the fact that xa is a projection. Non-commutativity of propo- 
sitions can be then trivially obtained if one notices that typical projectors of 
vectors do not commute (even in the 2D real plane the only projections that 
commute are those on parallel or perpendicular directions). 

So let now X be a vector and Pa a projector on direction A. Pa may be 
imagined as a matrix satisfying P^ — Pa- Let (-^^l^^) be a scalar product of 
two vectors. Take a unit X and define n{X) = {X\X) = 1, and 



{X\Pa\X) is a matrix element of Pa- If X^a-^^ ~ ^ (^^® identity matrix) 
then 



The set of projectors that are complete, i.e. that sum to I, represents a 
maximal collection of simultaneously "askable" questions. The associated 
probabilities sum to unity. However, even in the 2D plane there are infinitely 
many such sets. It follows that there may be infinitely many complete maxi- 
mal sets of simultaneously meaningful questions/propositions, but questions 
belonging to two different complete sets cannot be simultaneously considered 
(but can be asked one after another). 



The rule for joint and conditional probabilities becomes (if "first A then 
B" are "asked") 



p{A) = {PaX\PaX) = {X\Pa\X). 



(20) 




(21) 



A 



A 




{X\PaPbPa\X) = {X\EAnB\X). 



(22) 
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Probability in general depends on the order of questions, 

p{B r\A) = {X\PbPaPb\X) = {X\EBnA\X). (23) 

In symbolic notation this means that A ^ B occurs with different probability 
than B ^ A, a generic property of chemical or biological processes. Formula 



(22) represents probability of the answer "yes" to the question B, if it is 
asked in the context of a positive answer to an earlier question A. Note that 
EsnA and Eahb are projectors only if PaPb = PbPa- A general EsnA is 



the so called positive operator valued measure, POVM (Busch, Grabowski 



& Lahti, 1995). A POVM that is not a projector can be regarded as a 
non-Kolmogorovian analogue of a fuzzy-set membership function. 

The Kolmogorovian formula can be written as a commutative version of 
the same rule 

p{AnB) = ^^(x^X)^^^^^^^ = f^ixAXBXAX) 
= KxbXaXbX) = p{B n A). 

What we have sketched is just an example of a non-Kolmogorovian model. 
Another model, more directly applicable to population dynamics, is described 
in the next section. 

4. Density operator model of probability 

Let us employ Dirac's notation where complex column vectors X are 
denoted by |X). The scalar product 

n 

{Y\X) = J2y,X, = {Y„ . . . ,Yn) i (24) 
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can be regarded as a matrix multiplication of the 1 x n matrix 

{Y\ = {Y,,...Xn) = {\Y))^ (25) 



times the n x 1 matrix |X) (f denotes Hermitian conjugation). (24) regarded 
as a matrix product is, strictly speaking, not a number but a complex 1x1 
matrix (|y))^|X) = (y||X). Dirac's notation is based on identification of 
1x1 matrices with complex numbers, i.e. (y||X) = (Y\X). Now we can 
treat the formula p{A) = {X\Pa\X) as the product of three matrices: 1 x n 
{X\, n X n Pa, and n x 1 \X). The density operator formalism naturally 
occurs if one introduces the n x n matrix 



Px = \X){X\ 



fxA 



(Xi, . . . , Xn) 



( 



XiXn 



\Xri 



(26) 



satisfying 



p{A) 
Trpx 

{y\px\Y) 



(^I^aIx) 

Tr [PaPx] (definition) 
(X|X) = 1 (normalization) 
|(r|X)|2>0 (positivity) 
px (reality). 



(27) 
(28) 
(29) 
(30) 



Normalization, positivity and Hermiticity are characteristic of any combina- 
tion p = J2x ^xPx if Ax are probabilities, i.e. nonnegative real numbers 
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that sum to 1. Joint and conditional probabilities for a general p are defined 

by 

p{AnB) = Ti (PaPbPap) 

p{A) ^ V ^ 

piB\A) 

Hermiticity of p implies that its eigenvalues are real. Positivity means that 
these eigenvalues are nonnegative, and normalization guarantees that they 
sum to 1. If one skips normalization but keeps Hermiticity and positivity 
then for a projector P the number Tr (Pp) is nonnegative. If p{t) is a pos- 
itive Hermitian solution of some differential equation, then Tr (Pp(t)) is a 
nonnegative (kinetic) variable. 

Let us finally show that a single density operator p encodes in an ex- 



tremely efficient way a number of kinetic variables (Aerts & Czachor, 2006). 
p is an operator that acts in a linear (Hilbert) space H whose basis is given 
by a set {|^), (n|m) = 6nm} (there are infinitely many such bases). Matrix 
elements of p are in general complex 

Pnm = {n\p\m) = Xnm + iynm (32) 
Pnm ■^nm ^ynm •^mn ~l~ ^Vmn Pmn (^3) 

and thus Xmn •^nmi Umn llnrn) Unn 0. 

The diagonal elements p„„ = x„„ are themselves probabilities since 

Pnn = {n\p\n) = Tt PnP = Pn = Xnn (34) 

where P„ = \n){n\. Now, let 

= ^{\j) + \k)), (35) 
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m = ^{\j)-m), (36) 



Pjk^m{jk\,P;,^\jk'){jk'\. Then 



Xjk = Pjk - ^Pj - ^Pk, (37) 
Vjk = P'jk - \pj - \pk, (38) 



e 

Pjk = Pjk + ip'jk- -^{Pj+Pk): (39) 



where pjk = Tr PjkP, p'jk = T^t^ PjkP- foUows that a single p encodes three 
famihes of probabihties: {pn}, {Pjk}, and {pj^}- They are associated with 
three famihes of projectors: {Pn}, {Pjk}: and {Pjk}- Additional relations 
between the probabilities follow from 

Xjj = Pj = Pjj - Pj : Vjj = = p'jj - pj (40) 

and the resulting formula pj = p'jj = Pjj/'2. Let us note that {Pn} is complete, 
i.e. Yln=i^ = I, Yln=i^ Pn = 1- Order to understand completness 
properties of {Pjk} and {Pjk} we introduce, for j < k, two additional types 
of vectors and their associated projectors: 

= (41) 

= ^i\j)+m), (42) 

P^j^ = \jk^){jk^\, P'j^ = \jk'^){jk'^\. The completeness relations for {Pjk} 
and {Pjk} follow from the formula 

Pjk + P^k = Pj + Pk-Pjk + P'jk- (43) 

For any three Hcrmitian matrices satisfying [A, B] = iC one can prove 
the standard-deviation uncertainty relation AAAB > ||(C)| where AA — 
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■s/{A^) - (A)2, {A) = Tr(Ap), etc. If >1 = P = PMs a one-dimensional 
projector then (P) = p is a probability and one finds AP = -p). Two 

propositions Pi and P2 arc complementary if AP1AP2 > e > 0. Comple- 
mentarity means that reducing error APi we inevitably increase the one for 
AP2. 

In order to show that Pj and Pjk are complementary we compute 

[Pj.pA = \{\m-m\)^ic (44) 

and (C) = i/fcj, 

^Pj{l-Pj)\lpjkO--Pjk) >\\ykjV (45) 

The variable y^j measures complementarity of Pj and Pjk- 

The probabilities inherent in a single p have been so far defined with 
respect to a fixed basis \n). Being arbitrary, the basis could be replaced 
by any other orthonormal basis \n). Repeating the construction we would 
then arrive at new sets of projectors, P„ = ^^^d so on. They would 

lead to new families of probabilities, = Tr P„p, etc. There is no a priori 
rule that privileges one basis, so all these probabilities can be meaningful. 
What is important, all P„, P„m, P'nm^ may be complementary to all P„, 
Pnm, P'nm- What this practically means is that performing measurements 
of, say, the probability pi2, we inevitably influence possible future results of 
the remaining complementary probabilities. Measurement of pi2 creates a 
nontrivial conteoct for measurements of many (perhaps even all) probabilities 
Pn, Pnm, P'nrn^ as wcU as for p'^^, Pi, and P2. 

Finally, let us make the trivial remark that the standard Kolmogorovian 
model is a special case of the density operator formalism. The corresponding 
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p is then a combination 

n 

P-Y^P^Pj (46) 

where Pj belong to the same maximal set. This is equivalent to restricting 
p to diagonal matrices 



^ Pi ^ 



(47) 



\ Pn ) 

with (Kolmogorovian) probabilities on the diagonal. 

Nonlinear rate equations of a generalized Lotka-Volterra type can be for- 
mulated in terms of p{t). What is interesting, the formalism involving proba- 
bilities encoded by means of p{t) automatically allows us to formulate nonlin- 
ear rate equations in the so-called von Neumann, Liouville-von Neumann, or 
Lax forms. The latter property turns out to be essential for soliton structures 
and thus opens a possibility of solving very complicated coupled systems of 
rate equations by soliton techniques. In soliton and non-Kolmogorovian con- 
texts it is most appropriate to speak of Lax-von Neumann forms of rate 
equations: "Lax", since it stresses the sohton aspect, and "von Neumann" 
since the model of probability is based on density operators. 

4.1. Formal definition of contextuality 

Consider two random variables A and B such that the joint probability 

p{A = a (1 B = b) is defined for all values a and b of A and B, respectively. 
We assume that measurements of A are performed first, i.e. ^4 is a context for 
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B. Conditional and joint probabilities are defined (in both Kolmogorovian 
and non-Kolmogorovian frameworks) by the Bayes rule 

p{A = anB = b) =p{B = b\A = a)p{A = a). (48) 

W define the probability oi B = b in the context of A by 

p(5 = 6|A) = ^p(v4 = anfi = 6) (49) 

a 

The model (or problem) is non- contextual if 

p{B = b\Ai) = p{B = b\A2) = p{B = b) (50) 

for all random variables Ai and A2. In the density-matrix projector model 
the rule reads 



p{B = 6|A) = ^ Tr {pPA=aPB=bPA=a), 
a 

where ^^=a = T.b PB=h = I- If PA=aPB=b = PB=bPA=a then 

p{B = b\A) = Y,^^(pPA=aPB=bPA=a) 

a 

= 5^ Tr ipPA=aPB=b) = Tr (pP^^,) 

a 

= p{B = b). 



(51) 



(52) 



So, sensitivity to order of questions indicates contextuality. In (Aerts et 



al. , 2013) we show that contextuality in this sense is generic for nontrivial 



evolutionary games. Kolmogorovian models, with characteristic functions 
in the role of projectors, are non-contextual. Practical applications may 
require POVMs EA=anB=b yet more general than PA=aPB=bPA=a- This type 



of generalization is employed in (Aerts et al. , 2013) in order to reconstruct 



practically observed probabilities occurring in the RPS game played by Uta 
stanshuriana lizards. 
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5. Preliminaries on rate equations in Lax- von Neumann form 

Returning to the RPS game we can now solve the paradox. Namely, 
the dynamical aspect is localized in p{t) which is a matrix collecting all the 
possible propensities in all the possible contexts. In order to choose which 
context we need for S, say, we define two POVMs -Epns and -E^nSi so that 

Pnnsit) = TTp{t)Ejins, (53) 
ppnsit) = Ti p{t)Epns. (54) 

The dynamics of probabilities is defined through the dynamics of p{t). Con- 
textuality is present if J2l=0'^P=pns=s 7^ J2l=o '^R=rns=s- In Kolmogorovian 
probability the latter would be impossible since the characteristic functions 
Xpns and Xnns always do commute. It remains to define the dynamics of 
Pit). 

A Lax-von Neumann form of rate equations is 

p = -z[Hip),p]. (55) 

The dot denotes time derivative, [X, Y] = XY — YX, and H{p) is a linear 
Hermitian operator that depends on p. We assume that H{p) = H{p)'^ if 



p = p"^ . If pit) = r{pQ,t) is a concrete solution of (55) with initial condition 
p(0) = po, then 

H{p{t)) = H{r{po,t)) = h{po,t) (56) 

is a time-dependent operator satisfying, for this particular p(t) = r{pQ,t), 
the linear equation 

p{t) = -^[Kpo,t),pit)]. (57) 
18 



It is then known that there exists a unitary operator U{pQ,t) satisfying 

p{t) = U{po,t)p{0)U{po,t)l (58) 



The nonhnearity of (55) is reflected in the dependence of U{po,t) on the 



initial condition po- If U does not depend on po, i.e. is the same for all initial 



conditions then the dynamics t i— )■ p{t) is linear. The form of solution (58) 
is very important since it shows that p{t) and p(0) are related by unitary 
equivalence. Two unitarily equivalent Hermitian matrices have the same 
eigenvalues. In consequence, p{t) is positive whenever p(0) is positive. This, 
finally, guarantees that in order to guarantee positivity of Tr Pp(t) it is 
enough to start with a positive initial condition p(0). 

The central issue of the paper is the soliton, hence nonlinear dynamics of 
populations. However, since nonlinearities are not here exactly of the usual 
type, an interpretation in terms of abundances of populations vs. amounts of 
resources requires some new intuitions that are easier to develop on examples 
of linear evolutions. So let us begin with more detailed analysis of linear rate 
equations. 

5.1. Linear systems 

Let us now explicitly show a simple 2x2 example of standard rate equa- 
tions associated with their linear Lax-von Neumann form. In previous sec- 
tions we were not careful enough to distinguish between linear operators and 
their matrices, but now we need more precision. The same operator may 
be represented by different matrices, sice the notion of a matrix is basis de- 
pendent. So let h denote the matrix of the numbers hmn = {m\H\n), where 
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\n) are arbitrarily chosen orthonormal basis vectors. Similarly, let g be the 
matrix consisting of the numbers Qmn = {m\p\n). 
Now consider 

h = [ I. (69) 



Q21 Q22 



2 



(60) 



The Lax- von Neumann equation 

iQ = [h, q] (61) 
is equivalent to four rate equations, 

Pi = {ki2 - k'12) {Pi + P2) 

+2k[^pi2 - 2fci2p'i2> (62) 
P2 = {k'i2 - kl2) {Pi + P2) 

-2k[^Pi2 + 2A;i2Pi2, (63) 
P12 = (fci - ^2)^12 + ^- ^^ ^ - pi 

+ {-'^ + k',2) P2, (64) 

./ (ki-k2 \ fki-k2 \ 
P12 = ( — 2 — ^2 j ^ V — 2 j 

+ (A;2-fci)pi2. (65) 

It is clear that h is a matrix encoding kinetic constants of the dynamics. In 
order to identify appropriate types of interactions between the populations 
we have to know concrete values of the ks occurring in h. 
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So assume, for example, that 

ki = 1, k2 = 2ki, ki2 = 3/ci, k'12 = 4A;i, (66) 
and denote pi — A, p2 — B, pi2 — C, — Then 

A = -{A + B) + 8C -6D, (67) 

B = {A + B)-8C + 6D, (68) 

C = -l^+l^-D, (69) 

I) = l^-l^ + C. (70) 

We know by construction that A{t), B{t), C{t), D{t) are nonnegative if ^(0), 
-B(O), C(0), -D(O) are nonnegative, so this is a kinetic system, although not 
exactly of a chemical type. Adding the first two equations we find that 
A + B — Tr ^ is time independent (one of the conservation laws typical of 
Lax- von Neumann dynamics). 

Switching to another basis \h) we will obtain a different set of rate equa- 
tions. For example, the eigenvectors, H\n) = hn\n), make hmn = {m\H\n) 
diagonal, 

f hi \ 

h ^ { . (71) 

\ h2 J 

The corresponding probabilities 

Pn = Qnn = {n\p\n) = lY P„p = Tr {\n){n\p) (72) 



etc., satisfy 



P12 = k(pi2 — j, (73) 

P12 = -k[pi2 ^ — j, (74) 

Pi ^ p2^ 0, (75) 
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with k = hi — A general solution of these equations, 

Pi +P2 



+ (p'i2(0)- 



Pl2{t) 



-[pMo) 

Pi +P2 



Piit) 
P2{t) 



^Pi2(0) - 

Pi(0) = Pi, 

P2(0) = P2, 



P1+P2 

2 

Pi +P2 



P1+P2 

2 

Pi +P2 



COS kt 
sin H, 

sin kt 

cos H, 



(76) 



(77) 
(78) 
(79) 



with nonnegative initial condition at t = 0, remains nonnegative for all t 7^ 0. 
The sum of all the probabilities is not time independent since they do not 
belong to the same single maximal complete set. 



If we take the parameters (66), and denote A = pi, B = p2, C = pi2, 
D = p'12, then 

B = 0, (80) 

A + B^ 



A 
C 



D 

A + B 



A + B. 



2 

A + B 



(81) 

(82) 
(83) 



As we can see practically all elements of the kinetic system, including kinetic 
constants, have changed. 

The non-Kolmogorovity of the probability model allows for coexistence 



of (67)-(70) and (80)-(82) as representing context-dependent aspects of the 
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same dynamical system. One can show that there exists a time-independent 



terminology of information theory such a map is called a lossless communica- 
tion channel. These statements will hold true also for solutions of nonlinear 
equations discussed later on in this paper. 

5.2. Interpretation of the linear model 

All the dynamical variables are nonnegative so can be interpreted in terms 
of population abundances. The constant of motion (A + B)/2 = {A + B)/2 
defines a threshold that changes signs of derivatives of C{t) and D{t). Pop- 
ulation D{t) grows as long as the abundance C{t) is greater than {A + B)/2. 
When D{t) becomes greater than (^4 -|- -B)/2, the population represented by 
C{t) starts to decrease. At the moment the abundance C{t) becomes smaller 
than the threshold value {A + B)/2, the growth of D(t) stops, and a decay 
begins. 

Let us finally make the initial condition concrete. Assume A — B — 1/2, 



whose eigenvalues arc and 1. Since ^(0) is a positive operator (its eigen- 
values are nonnegative) then for any projector P one finds Tr {Pg{t)) > 0. 
In particular 



linear invertible transformation that maps A, B, C, D into A, B, C, D. In 



C(0) = 1, D{0) = 1/2, corresponding to 




(84) 



Tr(Pi2^(t)) = C{t) = - + - cos kt, 
Tr {P[,g{t)) = D{t) =^--^- sin kt. 



(85) 



(86) 
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C{t) and D{t) are complementary (belong to different complete sets) because 
[Pi2, -Pi'2] 7^ 0) s-^d thus do not have to sum to 1. They are shifted in phase 
with respect to each other in analogy to typical predator-prey abundances 
following from Lotka-Volterra models. 

6. Dynamics of subsystems in veirying environments 

Distinction between subsystems and their environments can be formalized 
in analogy to what one does in physics of open systems: State of a subsystem 
is influenced by the state of its environment, but not vice versa. Subsystems 
are coupled to environments in a non-symmetric way, a fact expressing the 
intuition that environments are "large" as compared to their inhabitants. If 
one cannot neglect the influence of a subsystem on its environment, then the 
whole "subsystem plus its environment" has to be considered as a single sys- 
tem, so that separation into two distinguished parts is no longer meaningful. 

Let us now consider the following coupled set of rate equations 

^1 = Hqi), (87) 

02 = F2{gi,g2), (88) 

^3 = -^3(^1,^2,^3), (89) 

; (90) 

Some of them may be representable in a Lax-von Neumann form, some of 
them perhaps not. The system described by ^1 plays a role of environment 
for the remaining subsystems. The one described by ^2 is the environment 
for Qs, g^, and so on. The collection of rate equations has to be solved in a 
hierarchical way. One begins with ^1 since the associated differential equation 
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is closed. Once one finds a given gi(t) = r(t), one switches to 

92 = F2{r,Q2). (91) 

Technically speaking, what one has to solve will be a set of coupled nonlinear 
rate equations with time dependent coefficients. At a first glance the problem 
looks, in its generality, hopelessly difficult. However, we will see that the 
power of soliton Lax- von Neumann equations may allow us to find explicit, 
exact, and highly nontrivial special solutions for the whole hierarchies of 
environments. 



7. Nonlinear systems 

Linear Lax-von Neumann equations lead to linear systems of rate equa- 
tions. Nonlinear rate equations of a generalized Lotka-Volterra type will 
occur if one takes less trivial H{g). For example, for H{g) = Ag + gA one 
finds 



ig 



[Ag + gA, g] = [A, g" 



(92) 



In the two-dimensional case, the nonlinearity in (92) is yet "too weak" since 



all nonlinear terms cancel out in the corresponding rate equations (this does 



not happen if A is a 3 x 3 matrix). (92) is a particular case of 



g = -i[HJ{g)] = -i[Hf{g),g] 



(93) 



where / is an arbitrary function and if is a linear operator. Given / one 
can find Hf{g), so this is indeed a Lax-von Neumann equation in the sense 



we have discussed above. The fact that (93) is a soliton rate system was 



estabhshed by Ustinov et al. (2001). 
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Yet another class of non-linear rate equations is obtained if one takes 



and replaces ( 71 ) by 



Hig) 



K 



k 
k 





2pi2 -Pi-P2 











The Lax-von Neumann equation 



(94) 



(95) 



ig = [H{g), g] 



(96) 



becomes equivalent to a system of coupled catalytic/auto-catalytic rate equa- 
tions, 



^, f Pi+P2\f , P1+P2 
P12 = 2ki^pi2 — )[Pi2- 



P12 



2k [pu 



Pi +P2V 



Pi = P2 = 0. 



(97) 

(98) 
(99) 



8. Soliton rate equations 

Soliton rate equations are not widely known but for details of their theory 
we refer the readers to the literature of the subject, beginning with the first 



work of Leble & Czachor (1998) where a soliton technique of integrating (92) 



was introduced. Generalization (93) of (92) is at the top of an infinite hier- 



archy of more complicated equations systematically cataloged by Cieslihski 
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Czachor & Ustinov (2003). Let us stress that H may be a differential operator 
and g could be infinitely dimensional. If one relaxes Hermiticity constraints 
then soliton Lax-von Neumann equations turn out to contain a large variety 
of integrable lattice systems. 

The term "soliton" is understood here in the general sense of "those equa- 
tions that are solvable by soliton methods". In technical terms this means 
that there exist Darboux-Backlund-covariant Lax pairs whose compatibility 
conditions are identical to the soliton system in question. 



Let us illustrate the latter statement by the soliton system (93). The Lax 
pair can be given here in various forms, but the following one is sufficent for 
our purposes, 

= {g-fiH)(p^, (100) 

i^t, = -fiQ)Vt,- (101) 
A* 

Here (time independent) complex numbers and y?^ is a matrix (it 

can be a vector, i.e. a 1-column matrix), yj^ can exist if 

= {ig- fiiH)i^^ + {g- fiH)-f{g)(p^ 



is equivalent to 



Since gf{g) = f{g)g, the condition reduces to 

{tg-[HJ{g)])^^ = tfiH^^ (102) 

If if = and ig~ [H, f{g)] = 0, which is our Lax-von Neumann system, then 
(f^ may exist. 
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The crucial step is given by the following theorem on Darhoux covariance 
of the Lax pair: Let iPfj.{t) be a 1-column matrix which is a solution of the 



Lax pair for some g{t) satisfying the compatibility condition (93). Then 



g[l] = g + {fi-fi)[P,H] (103) 
I + ^P),(I + ^P), (104) 

p = (105) 



is also a solution of (93) 



The theorem allows us to find a new solution g[l] given some known 
solution g. But how to find a, gl g cannot be too simple, say = I, or a more 
general g but commuting with H, since then ^[1] = g. Various tricks leading 
to nontrivial ^[1] have been nevertheless invented. For example, if f{g) = g^ 
one finds 

[H,g'] = [Hg + gH,g\. (106) 

A time independent g that anticommutes with H , gH = —Hg, leads to P 
that does not commute with H, so that g[l] ^ g. The procedure can be 
iterated: g = g\^ ^ g\\\ ^ g\2\ ^ . . .. In soliton terminology g\\\ is a 
1-soliton solution, since it is derived by means of a single Darboux-Backlund 
transformation, f)[0] — ?■ g\\\. 

I should be stressed that although dimensions of the matrices considered 
in our examples will be small, the method works in any dimension (even 
infinite). This is why Lax- von Neumann rate equations are naturally suited 
for problems involving multiple competition of a large number of species. If 
continuous time is replaced by a discrete time-step, the Lax-von Neumann 
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system turns into a kind of intransitive network whose vertices are defined 
by matrix elements of p. Intransitive networks liave recently been applied to 



the problem of biodiversity by AUesina & Levine (2011) 



So let us return to the problem of systems interacting with environments, 
and let 

Qi = F^{q^) (107) 

be any system of linear or nonlinear rate equations that describe the envi- 
ronment. In order to specify the form of F2 in 

Q2 = F2{gi, 02) (108) 

let us assume that in the absence of environment this is a soliton rate system 
involving (auto-)catalytic processes of order not higher than 2. Assuming 
that system-environment interaction also involves only elementary second- 
order processes we should finally obtain rate equations of order not higher 
than 3. A simple model possessing these characteristics is 

ig2 = ia + PTi{P,2gi))[A,il-s)g2 + sgl], (109) 

where a,P,s G M. Probabilities Pi^i2 = Tr(Pi2^i) are derived from the 
solution gi{t). A is an x n Hermitian matrix (ra = 3 is the lowest dimension 
for which quadratic terms do not cancel out in resulting rate equations). 
Parameter s allows us to control additional interactions between populations 
that were evolving independently of one another in the linear case discussed 
in Sec. 15.11 

We begin with the observation that if we know a solution w of 

iw = [A, {1 - s)w + sw% (110) 
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then 



g{t) = w{at + /3 /g rfrpi_i2(r)) 



;iii) 



solves (109). w{t) can be found by soliton techniques introduced by Leble & 



Czachor (1998). 



A relatively simple particular 1-soliton solution w = q[1] of (110) can be 
found if 

^ 
1 
2 

and we reduce the number of independent variables by imposing the con 



A 



\ 



(112) 



straints Wi2 = W23, wu = W33. Then (|Leble & Czachor[ 1998[) 

/ 

1 



wit) 



15 + 



\ 



;ii3) 



with 



at) 



{2 + 3i- V5t) VsTyl 



iult 



3(e27(t-to) + 1) 



;il4) 



;ii5) 



solves the Lax- von Neumann equation (110) (the readers may cross-check 



it directly in Wolfram Mathematica). The parameters are u = 1 



7 



15- 



and to £ is arbitrary. 



The associated set of rate equations for pa = P125 Vb 
Pd = Pi3 reads 



Pl3, PC = Pl2^ 



Pa 



-Pc- 



1-pi 1-pi 



Pi + 



30 



1 — Pi 1 — Pi 
+s I PiPa 7. — Pb H 7: — Pd 



Pb 



+s[pbPc -PaPd 

2pi - 2pD - 4:SpAPc 

(pi + l)(l-3p0 



(116) 
(117) 



PC 



+2s(-pi 

+2s (^{1 - pi)pa + (1 - Pi)pc + P2PD 
1 -pi 



+ PA + s [paPb + PcPd 
1 - Pi ^ 1 - Pi 



2 

+g(3 ^^ Pi-'^PiPA 2 



s PiPc 



1 -pi 



-Pd 



(118) 



Pd 



-2pi + 2pB + 2s {p2Pi - (1 - Pi)pa) 
+2s - P2PB + Pi)pc^ 
+2s(p^-p^). 



(119) 



For this concrete solution the diagonal elements are constants satisfying pi = 
P3) Pi -l- P2 -l- Ps = Trw = 1. Additional constants of motion are given 
by Tr (w") and Tr(Aw"'), for any natural n. An explicit solution of this 



complicated auto-catalytic system can be extracted from (113). Simplicity 



of the equivalent Lax- von Neumann equation (110) as compared to (116)- 



(119) is striking. 



Beginning again with the linear case s = we can interpret the black and 
green curves at Fig. 4 as resources for the consumers evolving according to 
the red and yellow curves (compare analogous plots for linear and nonlinear 



consumers given by Wilson & Abrams (2005)). Increasing s we gradually 
eliminate activity of the "red" consumer which approaches a steady state. 
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9. Periodicity of environment 



Seasonal changes can be modeled by periodic environment which couples 
to the two consumers and their resources described before. Periodicity of 
the environment, in its simplest version, can be modeled by a linear and 
conservative Lax- von Neumann dynamics of qi. In spite of high complication 
of the system of coupled nonlinear rate equations, let us stress again that we 
deal here with exact 1-soliton solutions. The reader may extract their explicit 



forms by means of (113), (111), and (85). 



Pb 



PC 



k[p'i,i2 
-Pc 



Figs. 7-11 show solutions of the coupled system 

Pl,12 
P'l,12 

Pa 



27' 
1 



12 



2 

1-pi 



■5 — ^ — Pi + 



+ S PiPA 



Pi 



-Pb 



Pi 



+s[pbPc - PaPd 



{a + (3pi, 



-PD 



12 



{a + /3pi,i2) 
+2s( -pi - 



2pi - 2pD - 4:spAPc 
■(pi + l)(l-3pi)' 



+2s (^{1 - Pi)pa + (1 - Pi)Pc + P2Pd 
1 -Pi 



^ +PA + S ypAPB + PcPd 
1 — Pl 1 — Pi 

+s{ 3 — - — Pl - 2pipA — Pb 



{a + /3pi,i2), 



I — Pl 
-slpipc H — Pd 
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Pd 



- 2pi + 2pB + 2s (^p2Pi - (1 - Pi)Pa^ 
+2s - P2PB + (1 - Pi)Pc^ 
+2s(p\-pI^ (a + /3pi,i2). 



Initial conditions for probabilities shown in Figs. 5-11 are identical in all the 
figures. We can say that the system is uncoupled from its environment if 
a — 1, /3 — 0. Parameters corresponding to a + /3/2 — lead to periodic 
dynamics of the system. Fig. 7 shows the effect of periodic "hibernation" of 
populations that can occur in seasons of small 

a + /?Pi,i2 ^ a + - + - coskt. (120) 

Figs. 7-8 involve strong coupling with environment, and two slightly differ- 
ent values of s. Parameter to is set to in all these plots, but we have 
the additional possibility of performing a controlled modification of initial 
conditions, if we modify to. Note that this would not be equivalent to just 
shifting the plots by to- All density matrices differing by the value of to be- 
long to the same "symplectic leaf of the dynamics, i.e. they posses the same 
eigenvalues. 

Fig. 12 shows sensitivity of the dynamics with respect to small changes of 
initial conditions. We again consider the case of no coupling to environment 
{a — 1, /3 — 0, s — Sc + 0.1). The initial conditions correspond to to — 30 
(solid curves) and to — 100 (dashed curves). The differences ai t — are 
smaller than thickness of the plots. 
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10. Non- normalized and decaying solutions 

It is clear that the standard context of population dynamics does not re- 
quire kinetic variables to be interpretable in terms of probabilities (positivity 
is sufficient). Moreover, typical rate equations involve mortality and birth 
rate terms that are absent in equations we have discussed so far. 

Let us briefly discuss these two issues. We begin with 

iw(t) = {l-s)[A,w(t)]+s[A,w(t)% (121) 
and define a new density operator 

v{t) = e'^^-'^^'w{t)e-'^^-'^^' (122) 



satisfying 



Now let 



iv{t) = s[A,v{tf]. (123) 



u{t) = x{t)v{y{t)) (124) 
where x{t) and y{t) are arbitrary differentiable functions. Then 

u{t) = i|y«W + |§(-OM,«(t)^]- (125) 

Finally, the unnormalized but positive matrix 

g{t) = e-^(^-^)^*xi(i)e^(i-^)^* 

= x(t)e^(^-^)^(^(*)-*V(|/(t))e-^(^-^)^(fW-*) (126) 

satisfies 

ig={l- s)[A, g] + -s[A, g'^] + i-g. (127) 

X X 
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If we are interested in rate equations with no explicit time dependence of 
parameters, we have to assume that 

x{t) = Xoe-^\ y{t) = yoe-^\ (128) 

leading to 

g = - s)[A, g] + t^[A, g^] - fxg. (129) 

Choosing, for example, Xq = —fiyo, we obtain the equation we have studied 
before, but with a term describing death (or birth) process, 

ig = {1 - s)[A, g] + s[A, g'^] - ing, (130) 

which is equivalent to 

Px = -fJ-Px + (131) 



X = A, B,C, D and the dots stand for the right-hand-sides of (116)-(119). 



The next three figures show solutions of ( 131 ) for three different couplings 
with environment. In all these figures we have the same initial conditions, 
k = 1, to = 0, fi = 0.1, s = 1, Xq = 1, yo = — Xo//i. In Fig. 13 the coupling 
parameters are a = 1, /3 = 0, meaning that the system is unaffected by its 
environment. In Fig. 14 a = 1 but /3 = —2. Finally, in Fig. 15 a = 1 and 
/3 = -1.9. 

As we can see, system that would exponentially decay in the absence 
of interaction with environment, may survive longer and even eliminate the 
decay if the coupling with environment is strong enough. 

It seems an appropriate place to mention the problem of paradox of the 



plankton (Hutchinson, 1961) or, more generally, the issue of biodiversity. 
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The phenomenon of species oscillations combined with nonlinear feedbacks 



(Huisman & Weissing, 1999) is a candidate for theoretical explanation of the 



failure of competitive exclusion principle (Gause, 1935). Here, periodically 



changing environment induces a similar phenomenon: The decay is either 
slowed by a kind of hibernation or even completely blocked by interaction 
with environment. 

So far we have included only the simplest form of dissipation, where all 
the right-hand sides of the equations are modified by the same linear term 
fip. One can further elaborate the theory by including the so called Lindblad 
terms of the form Ap = 2MpM^ - M^Mp - pM^M for some linear operator 
M. 

11. von Neumann form of generalized replicator equations 



Evolutionary game theory (Smith, 1982 Hofbauer & Sigmund, 1998) con 



centrates on evolution of probabilities of strategies. But any game involves 
two types of probabilities, existing at two essentially different levels. These 
second probabilities are implicitly present in matrix elements of payoff matri- 
ces. The latter probabilities do not couple to probabilities of strategies and 
play a role of parameters that determine the game. This is what happens at 
least in standard games, like poker. However, in games where players try to 
manipulate also the packs of cards, by adding aces hidden up their sleeves, 
the two levels get mixed with one another. Replicator equations apply to 
standard games, but aces-up-one 's-sleeves games must involve a generalized 
formalism. 



In (Aerts et al. , 2013) we show that the "pack of cards" probabilities 
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are in general non-Kolmogorovian. One can say that each entry of a payoff 
matrix involves probabilities associated with a separate probability space. 
An equation that mixes the two levels has to be based on a formalism that 
goes beyond Kolmogorovian probability. The Lax-von Neumann form turns 
out to be especially useful. 

The two types of probabilities enter the standard Kolmogorovian replica- 
tor equation, 

dx " " 

= X,[Y.-'^,X,-J2<^lmX,-m), (132) 

1=1 l,m=l 

N 

aim = "^bjpimj, /,m = l,...,n; j = l,...,N, 

n N 
^Xfc = 1; '^pimj = l; Xk>0,pimj>0, 
k=l j=l 

in different places. The non-Kolmogorovian ones are denoted by pim-j- The 
indices /, m in pim-j index probability spaces in a probability manifold, pimj 
are the "pack-of-cards" probabilities implicitly present in any game, but 
treated as fixed and independent of the probabilities of strategies Xk- The 
fact that the replicator equation possesses an equivalent Lax-von Neumann 



form was discovered by Gafiychuk & Prykarpatsky ( 2004 ) . Their form reads 

= [^i(Pi)>Pi], 
Pi = 

(1| = {^/xi,...,^/x:^), Xk=TrpiPl,k, 



n n 

,aiixi,. . . ,'^anixi ). 

1=1 1=1 



1 " 
-diag( J], 
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Here Pi^k = diag(0, . . . , 0, 1, 0, . . . 0), with 1 in the kth position. 

In order to generahze the Gafiychuk-Prykarpatsky equation to games with 
aces up one's sleeves one proceeds as follows. First denote puj = Pj, (1| = 
(a/P1) • • • ) a/Ptv), and p2 = It can be shown that there exists Eimj 



such that pim-j = p2Eim-y Elm-j = UimP2,jU}^ for some projectors P2 
Pj = Trp2-P2,j- Denote by Cg) the tensor product, and let p = pi ® p2, 
aim = Y.j bjh ® Ei-ra-j. Now, Xk = Trp(Pi,fc (g) I2), pim-j = Trp(Ii (g) Eim-j), 
aim = Tr pa/„. Denoting 

D{p) = -diag( J]aiiTrp(Pi,i®l2),...,^a„;Trp(Pi,i®l2)) ®l2, 

1=1 1=1 
H{p) = t[D{p),p], 

we reconstruct the standard replicator equation by taking the partial trace 
over the second subsystem from both sides of 

= [^(P),P], (133) 

under the constraints p = pi ® p2, dp2/dt = 0. The constraints imply 
[D{p),p] = [Di(pi), pi] ®p2, dp/dt = dpi/ dt® P2. Relaxing the constraints in 



(133), one generalizes (132) to games involving players with an ace up their 
sleeve, where correlations between and pim;j are no longer ignored. 

The question if replicator equations can be reformulated as a soliton sys- 
tem is an open one. 

12. Final remarks 



We have tried to show the potential inherent in soliton von Neumann 
equations. We have concentrated on the rate-equation aspect of the dynam- 
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ics. However, if matrices occurring in H are replaced by differential opera- 
tors one arrives at dynamical variables of a density type. The von Neumann 
equations then turn into a kind of reaction-diffusion systems, in general of 



an integro-differential type (Aerts et al. , 2003). The issues of dissipative 



systems, such as those preliminarily discussed in Section 10, require more 
detailed studies. Yet another possibility of including dissipation (based on 
a complex-time continuation of real-time soliton solutions) can be found in 



(Aerts & Czachor, 2007). The ace-up-one'e-sleeve generalization of replicator 
equation will be applied to a certain class of lizard-population evolutionary 



games in a forthcoming work. In the accompanying paper (Aerts et al. , 2013 ) 
we apply some of the ideas discussed here to a real life game played by Uta 
stanshuriana lizards. 
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Figure 1: Logical and probabilistic non-classicality of the rock-paper-scissors game, (a) 
There is no problem with three independent games. Results [S = Q,R = 1), {P = 0,S = 
1), (-R = 0, P = 1), occur with certainty, (b) Simultaneous game with three players is 
internally inconsistent. 
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PaPbX - 

Figure 2: Generic non-commutativity of "questions" represented by projections of vectors. 
The first "question" creates a context for the second one. The same "question" leads to 
context-dependent probabihties of answers. 




Figure 3: Linear model of population C (black) feeding on the resource D (red). 
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Figure 4: Dynamics (116)-(119) corresponding to s = 0, tg = 0. We observe two indepen- 



dent consumers pc and pu feeding on independent resources pa and pg . The environment 
acts trivially since a — 1, (5 = 0. For s = the dynamics is linear. 
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Figure 5: Dynamics (116)-(119l corresponding to s = 5 and the same initial condition as 



in Fig. 4. The environment still acts trivially since a = 1, j3 — 0. One of the consumers 
approaches steady state solution. 
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Figure 6: Dynamics corresponding to the critical value s — for which u) — 1— ^^^t^s^ = 
0. The same initial condition as in Fig. 4 and Fig. 5. The environment still acts trivially 
since a — 1, /3 — 0. The system approaches a steady state. 
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Figure 7: The same parameters as in Fig. 6, but now the systems are coupled to an 



oscillating environment (85): a ~ 0, (3 = 1, k = 1. We observe seasonal "hibernation" of 



all the four populations. The dotted curve is the oscillatory probability pi,i2(i)- 
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Figure 8: The same initial condition as before, but strong coupling to environment a = 5, 
^ = -9, fc 1. s = Sc - 0.025. 
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Figure 9: The same dynamics as in the previous figure, but for a longer time. We do not 
show the environment since its oscillation is too rapid in such a time scale. The next plot 
shows what happens if we increase s by 0.05... 
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Figure 12: The "butterfly effect". Sensitivity to initial conditions (no coupling to envi- 
ronment). Curves of the same color represent the same kinetic variable, but with slightly 
different values at < = 0. 
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Figure 13: Mortality rate is /i = 0.1. The remaining parameters are fc = l, ^o = 0, s = l, 
xq = 1, Ho = —xo/fi, a = 1, and /3 = 0. The latter means that the system is unaffected 
by its environment 
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Figure 14: The same as in the previous figure, but now /3 — —2, so that the environment 
is coupled to the system. The decay has been stopped due to periodicity ol environmental 
changes. The effect bears some resemblance to the paradox of the plankton. 
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Figure 15: The same parameters as in Fig. 14, but with /3 = —1.9, i.e. a slightly weaker 
coupling to environment. 
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